#1.筛选数据====================================================================
# 加载必要的包
library(tidyverse)
library(readr)
library(tidyr)
# 读取数据，注意文件编码
data_raw <- read_csv("SUOSUO.csv", locale = locale(encoding = "UTF-8"))  # 原始数据可能是GBK编码，根据实际情况调整

# 筛选梭梭，提取年份前四位
data_suoso <- data_raw %>%
  filter(植物名称 == "梭梭") %>%
  mutate(year = substr(测定日期, 1, 4)) %>%
  filter(year %in% c("2015", "2020", "2025")) %>%
  dplyr::select(height = 高度cm, crow = 冠幅cm, diam = 株丛径cm, year)

# 保存为TOTO.csv
write_csv(data_suoso, "TOTO.csv", na = "")


#补充数据
#=========================================================================
data_need <- read_csv("TOTO.csv")

# 定义填补函数
fill_na_diam <- function(df) {
  df_filled <- df
  na_idx <- which(is.na(df$diam))
  non_na <- df[!is.na(df$diam), ]
  
  for (i in na_idx) {
    h_target <- df$height[i]
    c_target <- df$crow[i]
    # 寻找高度和冠幅均在80%~120%范围内的观测
    candidates <- non_na %>%
      filter(height >= h_target * 0.8 & height <= h_target * 1.2,
             crow >= c_target * 0.8 & crow <= c_target * 1.2)
    if (nrow(candidates) > 0) {
      df_filled$diam[i] <- mean(candidates$diam, na.rm = TRUE)
    } else {
      # 若无，则用所有非缺失的均值
      df_filled$diam[i] <- mean(non_na$diam, na.rm = TRUE)
    }
  }
  return(df_filled)
}

data_filled <- fill_na_diam(data_need)

# 检查是否还有缺失
summary(data_filled)

# 保存
write_csv(data_filled, "trt_total.csv", na = "")


#龄级分类
#=====================================================================

library(readr)
library(dplyr)

# 读取数据
data <- read_csv("trt_total.csv")

# 先分别处理各年和整体
years <- c("2015", "2020", "2025")
data_list <- list()

for (yr in years) {
  sub <- data %>% filter(year == yr)
  # 计算年内最大值
  max_h <- max(sub$height)
  max_c <- max(sub$crow)
  max_d <- max(sub$diam)
  sub <- sub %>%
    mutate(importance = (height/max_h + crow/max_c + diam/max_d)/3 * 100,
           year_group = yr)
  
  # 按株丛径最大值最小值等间距划分10级（替换原百分位数逻辑）
  min_d <- min(sub$diam, na.rm = TRUE)  # 年内diam最小值
  max_d_year <- max(sub$diam, na.rm = TRUE)  # 年内diam最大值
  sub <- sub %>%
    mutate(diam_quant = cut(diam, 
                            breaks = seq(min_d, max_d_year, length.out = 11),  # 等间距生成11个分割点（10级）
                            include.lowest = TRUE, labels = FALSE)) %>%
    rename(diam_class_year = diam_quant)
  
  # 按重要指数最大值最小值等间距划分（替换原百分位数逻辑）
  min_imp <- min(sub$importance, na.rm = TRUE)  # 年内importance最小值
  max_imp <- max(sub$importance, na.rm = TRUE)  # 年内importance最大值
  sub <- sub %>%
    mutate(imp_quant = cut(importance,
                           breaks = seq(min_imp, max_imp, length.out = 11),  # 等间距生成11个分割点（10级）
                           include.lowest = TRUE, labels = FALSE)) %>%
    rename(imp_class_year = imp_quant)
  
  data_list[[yr]] <- sub
}

# 合并三年数据
total <- bind_rows(data_list)
# 计算整体最大值
max_h_total <- max(total$height)
max_c_total <- max(total$crow)
max_d_total <- max(total$diam)
total <- total %>%
  mutate(importance_total = (height/max_h_total + crow/max_c_total + diam/max_d_total)/3 * 100)

# 整体按株丛径最大值最小值等间距划分（替换原百分位数逻辑）
total <- total %>%
  group_by(year) %>%  
  ungroup() %>%
  mutate(diam_quant_total = cut(diam,
                                breaks = seq(min(total$diam, na.rm = TRUE), 
                                             max(total$diam, na.rm = TRUE), 
                                             length.out = 11),  # 整体等间距10级
                                include.lowest = TRUE, labels = FALSE))

# 整体按重要指数最大值最小值等间距划分（替换原百分位数逻辑）
total <- total %>%
  mutate(imp_quant_total = cut(importance_total,
                               breaks = seq(min(total$importance_total, na.rm = TRUE), 
                                            max(total$importance_total, na.rm = TRUE), 
                                            length.out = 11),  # 整体等间距10级
                               include.lowest = TRUE, labels = FALSE))

# 重新整理列顺序
data_final <- total %>%
  dplyr::select(year, height, crow, diam, diam_class_year, imp_class_year, 
         diam_quant_total, imp_quant_total, importance, importance_total)

# 写入数据
write_csv(data_final, "data.csv", na = "")


#4.计算生命表
#============================================================================
# 生命表 + 存活/死亡率/消失率曲线（简化版）
# 核心功能：生命表计算、Deevey类型判断、SCI标准图表生成

library(ggplot2)
library(gridExtra)
library(readr)
library(dplyr)
library(RColorBrewer)
library(grid)

# 读取数据
data <- read_csv("data.csv", show_col_types = FALSE)

# 1. 生命表构建函数
build_life_table <- function(df, class_col, name) {
  df <- df %>% filter(!is.na(!!sym(class_col)))
  counts <- df %>% count(!!sym(class_col)) %>% arrange(!!sym(class_col))
  
  nx   <- counts$n
  lx   <- nx / sum(nx) * 1000  # 标准化为1000
  dx   <- c(lx[1:(length(lx)-1)] - lx[2:length(lx)], NA)
  qx   <- dx / lx
  kx   <- c(log(lx[1:(length(lx)-1)]) - log(lx[2:length(lx)]), NA)
  
  data.frame(
    age_group = counts[[1]],
    nx = nx, lx = lx, dx = dx, qx = qx, kx = kx,
    method = name
  )
}

# 2. 生成生命表数据
# Diameter
lt_diam_total <- build_life_table(data, "diam_class_year", "total")
lt_diam_year <- data %>%
  group_by(year) %>%
  do(build_life_table(., "diam_class_year", as.character(.$year[1]))) %>%
  ungroup()

# Importance
lt_imp_total <- build_life_table(data, "imp_class_year", "total")
lt_imp_year <- data %>%
  group_by(year) %>%
  do(build_life_table(., "imp_class_year", as.character(.$year[1]))) %>%
  ungroup()

# 数据清洗
life_diam <- bind_rows(lt_diam_total, lt_diam_year) %>%
  filter(!is.na(lx), !is.na(qx), !is.na(kx)) %>%
  mutate(method = factor(method, levels = c("total","2015","2020","2025")))

life_imp <- bind_rows(lt_imp_total, lt_imp_year) %>%
  filter(!is.na(lx), !is.na(qx), !is.na(kx)) %>%
  mutate(method = factor(method, levels = c("total","2015","2020","2025")))

# 保存生命表
life_all <- bind_rows(
  life_diam %>% mutate(type="diam"),
  life_imp  %>% mutate(type="imp")
)
write_csv(life_all, "FIG01.csv")

# 3. Deevey类型判断函数（简化版）
fit_survival_models <- function(plot_data) {
  x <- plot_data$age_group
  y <- plot_data$lx
  
  # 基础检查
  if(length(x) < 3) {
    return(list(code = "?", type = "Unknown", drop = 0, r2_exp = 0, r2_lin = 0, r2_power = 0))
  }
  
  # 计算幼体下降率
  juvenile_x <- sort(unique(x))[1:min(3, length(unique(x)))]
  juvenile_data <- plot_data %>% filter(age_group %in% juvenile_x) %>% arrange(age_group)
  
  drop_rate <- if(nrow(juvenile_data) >=2 && juvenile_data$lx[1] > 0) {
    (juvenile_data$lx[1] - juvenile_data$lx[nrow(juvenile_data)])/juvenile_data$lx[1]
  } else 0
  
  # 模型拟合
  exp_mod <- tryCatch(nls(y ~ a*exp(b*x), start=list(a=max(y)*1.2, b=-0.05),
                          control=nls.control(maxiter=2000), algorithm="port",
                          lower=c(0,-1), upper=c(max(y)*2,0)), error=function(e) NULL)
  
  lin_mod <- lm(y ~ x)
  
  pow_mod <- tryCatch(nls(y ~ a*x^b, start=list(a=max(y)*5, b=-0.5),
                          control=nls.control(maxiter=2000), algorithm="port",
                          lower=c(0,-5), upper=c(max(y)*10,0)), error=function(e) NULL)
  
  # 计算R²
  r2_exp <- if(!is.null(exp_mod)) {
    pred <- pmax(predict(exp_mod), 0)
    1 - sum((y-pred)^2)/sum((y-mean(y))^2)
  } else 0
  
  r2_lin <- summary(lin_mod)$r.squared
  
  r2_power <- if(!is.null(pow_mod)) {
    pred <- pmax(predict(pow_mod), 0)
    1 - sum((y-pred)^2)/sum((y-mean(y))^2)
  } else 0
  
  # 判断类型
  max_r2 <- max(r2_exp, r2_lin, r2_power)
  is_convex <- (r2_exp == max_r2) | (drop_rate < 0.3)
  is_concave <- (r2_power >= 0.6) | (drop_rate >= 0.3)
  is_linear <- (r2_lin - max(r2_exp, r2_power) > 0.1) & !is_convex & !is_concave
  
  if(is_convex) {
    code <- "Ⅰ"; type <- "Convex (Ⅰ type)"
  } else if(is_concave) {
    code <- "Ⅲ"; type <- "Concave (Ⅲ type)"
  } else if(is_linear) {
    code <- "Ⅱ"; type <- "Linear (Ⅱ type)"
  } else {
    code <- "?"; type <- "Unknown"
  }
  
  return(list(code = code, type = type, drop = round(drop_rate,3),
              r2_exp = round(r2_exp,3), r2_lin = round(r2_lin,3), r2_power = round(r2_power,3)))
}

# 4. 拟合Deevey类型
diam_res <- fit_survival_models(filter(life_diam, method=="total"))
imp_res  <- fit_survival_models(filter(life_imp, method=="total"))

# 5. SCI主题（简化版）
theme_sci <- function(show_legend=FALSE) {
  theme_bw() +
    theme(
      text = element_text(colour="black", size=10),
      plot.title = element_text(size=12, face="bold", hjust=0.5, margin=margin(b=5)),
      axis.title = element_text(size=11, face="bold", margin=margin(t=5)),
      axis.text = element_text(size=10, colour="black"),
      axis.line = element_line(linewidth=0.5, colour="black"),
      axis.ticks = element_line(linewidth=0.5, colour="black"),
      axis.ticks.length = unit(3, "mm"),
      panel.background = element_blank(),
      panel.grid = element_blank(),
      panel.border = element_blank(),
      legend.title = element_text(size=10, face="bold"),
      legend.text = element_text(size=9),
      legend.key = element_rect(fill="white", colour=NA),
      legend.background = element_rect(fill="white", colour=NA),
      legend.margin = margin(0,0,0,0),
      plot.margin = unit(c(0.8,0.8,0.8,0.8), "cm"),
      legend.position = if(show_legend) "right" else "none"
    )
}

# 配色
cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3")
names(cols) <- c("total","2015","2020","2025")

# 6. 图表绘制函数（合并重复逻辑）
plot_life_table <- function(data, type_name, deevey_res) {
  # 存活曲线
  p1 <- ggplot(data, aes(x=age_group, y=lx, group=method, color=method)) +
    geom_line(linewidth=1.2, alpha=0.9) +
    geom_point(size=2.5, shape=16) +
    scale_color_manual(values=cols) +
    scale_x_continuous(breaks=seq(1, max(data$age_group),1), expand=c(0.02,0)) +
    scale_y_continuous(expand=c(0.02,0)) +
    labs(x="Age class", y="Standardized survivorship (lx)",
         title=paste0(type_name, " – Survivorship curve (", deevey_res$type, ")")) +
    theme_sci()
  
  # 死亡率曲线
  p2 <- ggplot(data, aes(x=age_group, y=qx, group=method, color=method)) +
    geom_line(linewidth=1.2, alpha=0.9) +
    geom_point(size=2.5, shape=16) +
    scale_color_manual(values=cols) +
    scale_x_continuous(breaks=seq(1, max(data$age_group),1), expand=c(0.02,0)) +
    scale_y_continuous(expand=c(0.02,0)) +
    labs(x="Age class", y="Mortality rate (qx)",
         title=paste0(type_name, " – Mortality curve")) +
    theme_sci()
  
  # 消失率曲线（显示图例）
  p3 <- ggplot(data, aes(x=age_group, y=kx, group=method, color=method)) +
    geom_line(linewidth=1.2, alpha=0.9) +
    geom_point(size=2.5, shape=16) +
    scale_color_manual(values=cols, name="Time period") +
    scale_x_continuous(breaks=seq(1, max(data$age_group),1), expand=c(0.02,0)) +
    scale_y_continuous(expand=c(0.02,0)) +
    labs(x="Age class", y="Killing power (kx)",
         title=paste0(type_name, " – Killing power curve")) +
    theme_sci(show_legend=TRUE)
  
  # 总标题
  main_title <- textGrob(paste0(type_name, "-based Life Table Analysis"),
                         gp=gpar(fontsize=14, fontface="bold"), y=1, vjust=1)
  
  # 组合图表
  p_combined <- grid.arrange(p1, p2, p3, ncol=3, widths=c(1,1,1.2), top=main_title)
  
  # 保存
  ggsave(paste0("FIG02_", tolower(type_name), "_SCI.png"), p_combined,
         width=18, height=6, dpi=600, units="in", bg="white")
  ggsave(paste0("FIG02_", tolower(type_name), "_SCI.pdf"), p_combined,
         width=18, height=6, device=cairo_pdf)
  
  return(p_combined)
}

# 7. 生成图表
plot_life_table(life_diam, "Diameter", diam_res)
plot_life_table(life_imp, "Importance", imp_res)

# 8. 输出结果
cat("=== Deevey Type Fitting Results ===\n")
cat("★ Diameter: ", diam_res$type, " (幼体下降率: ", diam_res$drop, ")\n",
    "  R²: 指数(Ⅰ)=", diam_res$r2_exp, " 线性(Ⅱ)=", diam_res$r2_lin, " 幂指数(Ⅲ)=", diam_res$r2_power, "\n\n",
    "★ Importance: ", imp_res$type, " (幼体下降率: ", imp_res$drop, ")\n",
    "  R²: 指数(Ⅰ)=", imp_res$r2_exp, " 线性(Ⅱ)=", imp_res$r2_lin, " 幂指数(Ⅲ)=", imp_res$r2_power, "\n\n",
    "=== 文件输出 ===\n",
    "✓ 生命表: FIG01.csv\n",
    "✓ 图表: FIG02_diam_SCI.png/pdf, FIG02_imp_SCI.png/pdf\n", sep="")

#计算多样性指数
#====================================================================
library(vegan)  # 用于多样性指数

data <- read_csv("data.csv")

# 定义计算多样性的函数
calc_diversity <- function(df, class_col, group_name) {
  counts <- df %>% count(!!sym(class_col)) %>% pull(n)
  if (length(counts) == 0) return(NULL)
  p <- counts / sum(counts)
  H <- -sum(p * log(p))
  D <- 1 - sum(p^2)
  J <- H / log(length(p))
  return(data.frame(group = group_name, method = class_col,
                    Shannon = H, Simpson = D, Pielou = J))
}

# 定义所有要计算的组合
combinations <- list(
  list(df = data, class = "diam_class_year", name = "2015_diam"),
  list(df = filter(data, year == "2015"), class = "diam_class_year", name = "2015_diam"),
  list(df = filter(data, year == "2020"), class = "diam_class_year", name = "2020_diam"),
  list(df = filter(data, year == "2025"), class = "diam_class_year", name = "2025_diam"),
  list(df = data, class = "imp_class_year", name = "2015_imp"),
  list(df = filter(data, year == "2015"), class = "imp_class_year", name = "2015_imp"),
  list(df = filter(data, year == "2020"), class = "imp_class_year", name = "2020_imp"),
  list(df = filter(data, year == "2025"), class = "imp_class_year", name = "2025_imp"),
  list(df = data, class = "diam_quant_total", name = "total_diam"),
  list(df = data, class = "imp_quant_total", name = "total_imp")
)

# 注意：上面有些重复，需要修正。实际上我们需要：
# 对于每个年份，分别用diam_class_year和imp_class_year；对于整体，用diam_class_total和imp_class_total。
# 所以总共：2015_diam, 2015_imp, 2020_diam, 2020_imp, 2025_diam, 2025_imp, total_diam, total_imp
# 重新整理
div_results <- bind_rows(
  calc_diversity(filter(data, year == "2015"), "diam_class_year", "2015_diam"),
  calc_diversity(filter(data, year == "2015"), "imp_class_year", "2015_imp"),
  calc_diversity(filter(data, year == "2020"), "diam_class_year", "2020_diam"),
  calc_diversity(filter(data, year == "2020"), "imp_class_year", "2020_imp"),
  calc_diversity(filter(data, year == "2025"), "diam_class_year", "2025_diam"),
  calc_diversity(filter(data, year == "2025"), "imp_class_year", "2025_imp"),
  calc_diversity(data, "diam_quant_total", "total_diam"),
  calc_diversity(data, "imp_quant_total", "total_imp")
)

# 绘制三个指数的图
p_shannon <- ggplot(div_results, aes(x = group, y = Shannon, fill = method)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "", y = "Shannon-Wiener index", title = "Shannon diversity") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_simpson <- ggplot(div_results, aes(x = group, y = Simpson, fill = method)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "", y = "Simpson index", title = "Simpson diversity") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_pielou <- ggplot(div_results, aes(x = group, y = Pielou, fill = method)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "", y = "Pielou evenness", title = "Pielou evenness") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 合并三个图，或者分别保存
library(cowplot)
p_combined_div <- plot_grid(p_shannon, p_simpson, p_pielou, ncol = 3, labels = "AUTO")
ggsave("FIG03.png", p_combined_div, width = 15, height = 5, dpi = 300)
ggsave("fig03.pdf", p_combined_div, width = 15, height = 5)
write_csv(div_results , "div_results.csv")


#计算生态位宽度和重叠指数计算
#========================================================================
# 加载必要的包
library(tidyverse)
library(reshape2)
library(gridExtra)

# 读取数据
df <- read.csv("data.csv")

# ================== 1. 定义函数 ==================
# 计算Levins宽度和标准化宽度
levins <- function(counts, n_classes) {
  p <- counts / sum(counts)
  B <- 1 / sum(p^2)
  Ba <- (B - 1) / (n_classes - 1)
  return(c(B = B, Ba = Ba))
}

# 计算Morisita重叠指数
morisita <- function(p, q) {
  2 * sum(p * q) / (sum(p^2) + sum(q^2))
}

# 获取所有可能的龄级（资源状态）
diam_classes <- sort(union(unique(df$diam_class_year), unique(df$diam_quant_total)))
imp_classes <- sort(union(unique(df$imp_class_year), unique(df$imp_quant_total)))

# 定义分组，使用英文避免字体问题
groups <- c("2015", "2020", "2025", "Total")

# ================== 2. 直径分析 ==================
# 构建数据列表
diam_data <- list(
  "2015" = df %>% filter(year == 2015) %>% pull(diam_class_year),
  "2020" = df %>% filter(year == 2020) %>% pull(diam_class_year),
  "2025" = df %>% filter(year == 2025) %>% pull(diam_class_year),
  "Total" = df %>% pull(diam_quant_total)
)

# 计算各组的频数和比例
diam_counts <- lapply(diam_data, function(x) table(factor(x, levels = diam_classes)))
diam_props <- lapply(diam_counts, function(x) x / sum(x))

# 计算宽度指数
diam_width <- sapply(diam_counts, function(cnt) levins(cnt, length(diam_classes))) %>%
  t() %>% as.data.frame() %>%
  mutate(Group = groups, .before = 1)

# ================== 3. 重要性分析 ==================
imp_data <- list(
  "2015" = df %>% filter(year == 2015) %>% pull(imp_class_year),
  "2020" = df %>% filter(year == 2020) %>% pull(imp_class_year),
  "2025" = df %>% filter(year == 2025) %>% pull(imp_class_year),
  "Total" = df %>% pull(imp_quant_total)
)

imp_counts <- lapply(imp_data, function(x) table(factor(x, levels = imp_classes)))
imp_props <- lapply(imp_counts, function(x) x / sum(x))

imp_width <- sapply(imp_counts, function(cnt) levins(cnt, length(imp_classes))) %>%
  t() %>% as.data.frame() %>%
  mutate(Group = groups, .before = 1)

# ================== 4. 绘制柱形图 ==================
# 直径宽度柱形图
p_diam <- diam_width %>%
  pivot_longer(cols = c(B, Ba), names_to = "Index", values_to = "Value") %>%
  ggplot(aes(x = Group, y = Value, fill = Index)) +
  geom_col(position = position_dodge(0.9)) +
  labs(title = "Levins Niche Width (Diameter Class)",
       y = "Width", x = "Group") +
  theme_minimal() +
  scale_fill_manual(values = c("B" = "steelblue", "Ba" = "coral"))

ggsave("FIG04-diam.png", p_diam, width = 6, height = 4)
ggsave("FIG04-diam.pdf", p_diam, width = 6, height = 4)

# 重要性宽度柱形图
p_imp <- imp_width %>%
  pivot_longer(cols = c(B, Ba), names_to = "Index", values_to = "Value") %>%
  ggplot(aes(x = Group, y = Value, fill = Index)) +
  geom_col(position = position_dodge(0.9)) +
  labs(title = "Levins Niche Width (Importance Class)",
       y = "Width", x = "Group") +
  theme_minimal() +
  scale_fill_manual(values = c("B" = "steelblue", "Ba" = "coral"))

ggsave("FIG04-imp.png", p_imp, width = 6, height = 4)
ggsave("FIG04-imp.pdf", p_imp, width = 6, height = 4)

# ================== 5. 计算Morisita重叠指数 ==================
# 直径重叠矩阵
n_groups <- length(groups)
diam_overlap <- matrix(1, n_groups, n_groups, dimnames = list(groups, groups))

for (i in 1:(n_groups-1)) {
  for (j in (i+1):n_groups) {
    p_i <- diam_props[[i]]
    p_j <- diam_props[[j]]
    overlap <- morisita(p_i, p_j)
    diam_overlap[i, j] <- overlap
    diam_overlap[j, i] <- overlap
  }
}

# 重要性重叠矩阵
imp_overlap <- matrix(1, n_groups, n_groups, dimnames = list(groups, groups))

for (i in 1:(n_groups-1)) {
  for (j in (i+1):n_groups) {
    p_i <- imp_props[[i]]
    p_j <- imp_props[[j]]
    overlap <- morisita(p_i, p_j)
    imp_overlap[i, j] <- overlap
    imp_overlap[j, i] <- overlap
  }
}

# ================== 6. 绘制热图 ==================
# 直径热图
melt_diam <- melt(diam_overlap)
names(melt_diam) <- c("Group1", "Group2", "Overlap")
p_heat_diam <- ggplot(melt_diam, aes(x = Group1, y = Group2, fill = Overlap)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "red", limits = c(0, 1)) +
  labs(title = "Morisita Overlap (Diameter Class)", x = "", y = "") +
  theme_minimal() +
  coord_fixed()

ggsave("FIG04-diam-heatmap.png", p_heat_diam, width = 5, height = 4)
ggsave("FIG04-diam-heatmap.pdf", p_heat_diam, width = 5, height = 4)

# 重要性热图
melt_imp <- melt(imp_overlap)
names(melt_imp) <- c("Group1", "Group2", "Overlap")
p_heat_imp <- ggplot(melt_imp, aes(x = Group1, y = Group2, fill = Overlap)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "red", limits = c(0, 1)) +
  labs(title = "Morisita Overlap (Importance Class)", x = "", y = "") +
  theme_minimal() +
  coord_fixed()

ggsave("FIG04-imp-heatmap.png", p_heat_imp, width = 5, height = 4)
ggsave("FIG04-imp-heatmap.pdf", p_heat_imp, width = 5, height = 4)

# ================== 7. 保存计算结果为ECWR.csv（优化版，无警告） ==================
# 先写入直径宽度（包含列名）
write.table(diam_width, file = "ECWR.csv", sep = ",", row.names = FALSE,
            col.names = TRUE, quote = FALSE)

# 追加直径重叠矩阵（手动处理列名，避免警告）
cat("\n# Diameter Overlap Matrix\n", file = "ECWR.csv", append = TRUE)
# 写入列名行：第一列为空（对应行名），后面为分组名
cat(paste0(",", paste(colnames(diam_overlap), collapse = ","), "\n"), 
    file = "ECWR.csv", append = TRUE)
# 写入矩阵数据（不包含列名）
write.table(diam_overlap, file = "ECWR.csv", sep = ",", row.names = TRUE,
            col.names = FALSE, quote = FALSE, append = TRUE)

# 追加重要性宽度
cat("\n# Importance Width\n", file = "ECWR.csv", append = TRUE)
write.table(imp_width, file = "ECWR.csv", sep = ",", row.names = FALSE,
            col.names = TRUE, quote = FALSE, append = TRUE)

# 追加重要性重叠矩阵
cat("\n# Importance Overlap Matrix\n", file = "ECWR.csv", append = TRUE)
cat(paste0(",", paste(colnames(imp_overlap), collapse = ","), "\n"), 
    file = "ECWR.csv", append = TRUE)
write.table(imp_overlap, file = "ECWR.csv", sep = ",", row.names = TRUE,
            col.names = FALSE, quote = FALSE, append = TRUE)

cat("所有计算和绘图完成！结果已保存，且CSV文件已优化，无警告。\n")

#计算稳定性指数
#============================================================================
data <- read_csv("data.csv")

# 用重要指数（整体）划分的龄级 imp_class_total 来计算
stability <- data %>%
  group_by(imp_quant_total) %>%
  summarise(
    mean_imp = mean(importance_total),
    sd_imp = sd(importance_total),
    stab = mean_imp / sd_imp
  ) %>%
  ungroup()

# 计算整体稳定指数（所有龄级稳定指数的均值/标准差）
overall_stab <- mean(stability$stab, na.rm = TRUE) / sd(stability$stab, na.rm = TRUE)

# 对于各年份，也需要计算
stability_year <- data %>%
  group_by(year, imp_quant_total) %>%
  summarise(
    mean_imp = mean(importance_total),
    sd_imp = sd(importance_total),
    stab = mean_imp / sd_imp,
    .groups = "drop"
  ) %>%
  group_by(year) %>%
  summarise(
    overall_stab = mean(stab, na.rm = TRUE) / sd(stab, na.rm = TRUE)
  ) %>%
  ungroup()

# 合并整体和年份
stab_all <- bind_rows(
  stability_year,
  data.frame(year = "total", overall_stab = overall_stab)
)

# 绘图
p_stab <- ggplot(stab_all, aes(x = year, y = overall_stab)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(x = "", y = "Stability index", title = "Population stability") +
  theme_minimal()
ggsave("FIG06.png", p_stab, width = 8, height = 6, dpi = 300)
ggsave("fig06.pdf", p_stab, width = 8, height = 6)







































